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ABSTRACT 

The asymmetric emission of gravitational waves produced during the coalescence of a massive black 
hole (MBH) binary imparts a velocity "kick" to the system that can displace the hole from the center 
of its host. Here we study the trajectories and observability of MBHs recoiling in three (one major, 
two minor) gas-rich galaxy merger remnants that were previously simulated at high resolution, and 
in which the pairing of the MBHs had been shown to be successful. We run new simulations of 
MBHs recoiling in the major merger remnant with Mach numbers in the range 1 < A4 < 6, and use 
simulation data to construct a semi-analytical model for the orbital evolution of MBHs in gas-rich 
systems. We show that: 1) in major merger remnants the energy deposited by the moving hole into 
the rotationally supported, turbulent medium makes a negligible contribution to the thermodynamics 
of the gas. This contribution becomes significant in minor merger remnants, potentially allowing for 
an electromagnetic signature of MBH recoil; 2) in major merger remnants, the combination of both 
deeper central potential well and drag from high-density gas confines even MBHs with kick velocities 
as high as 1,200 kms -1 within 1 kpc from the host's center; 3) kinematically offset nuclei may be 
observable for timescales of a few Myr in major merger remnants in the case of recoil velocities in 
the range 700 — 1,000 kms -1 ; 4) in minor mergers remnants the effect of gas drag is weaker, and 
MBHs with recoil speeds in the range 300 — 600 kms -1 will wander through the host halo for longer 
timescales. When accounting for the probability distribution of kick velocities, however, we find that 
the likelihood of observing recoiling MBHs in gas-rich galaxy mergers is very low, typically below 

io -5 -io -6 . 

Subject headings: black hole physics - galaxies: halos - kinematics and dynamics - methods: numerical 



1. INTRODUCTION 



The pairing of massive black holes (MBHs) is a natu- 
ral consequence of galaxy mergers in hierarch i cal struc 
ture format ion scenarios ( [Begelman et al.J1980 



Volonteri 



ical relativity (Pretorius 2005 [Campanelli et al. 2006|) 
have allowed several gro ups (Baker et al. 2006; Her- 
rmann et al. 2007 Gonzalez et al.||2007[ ) to evolve bi- 
nary black holes from inspiral through coalescence and 
estimate the resultin g recoil velocities as a function of 



EFld1f2003|). To date, o ver 30 dual actiVe galactic mi- estimate the resulting recoil velocities as a function oi 
clei (JtUNs) have been discovered in merger remnants at ™ atl ^ and spin ( [Baker et al |2008| Campanelli et al. 



redshifts 0.34 < z < 0.82 in the DEEP Survey (Comer- 



ford et al. 112009), over 160 have be reporte d from st udies 



of double-peaked [PHI] profiles (e.g. |Zhou et al.||2004l 



Gerke et al. 



2007 Liu et al. 2010 Smith et al. 2010) 



and n early 40 from the Swift BAT survey ( |Koss et al.| 



2010 ). Furthermore, Chandra observations of the nucleus 
of 1NGC6240 have unveiled a MBH binary with se para- 



tion < 1 kpc (jKomossa et al.|2003||Max et al |2007| ). The 

coalescence of MBH binaries may become the dominant 
source of gravitational waves at mHz frequencies (e.g. , 
Haeh nelt|[l994l |Wyithe & Loeb 2003; Sesana et al] 2004 



2005p , a sensitivity range probed by the Laser Inter fer- 
ometer Space Antenna (LISA). 

Asymmetries in the configuration of the coalescing bi- 
nary cause the beaming of gravitational wave radiation 
in a preferred direction. To conserve linear momentum, 
the remnant h ole recoils in the direction opposite to the 
radiative flux ( l Peres|| 19621 |Bekenstein||1973| |Fitchett fc| 
|Detweiler||1984p . The recoil or "kick" velocity, Vkick, de- 
pends on the binary mass ratio and on the magnitude 
and direction of their spins, and does not depend on the 
total mass of the binary. Recent advances in numer- 
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2007a). These studies show kick velocities in the range 
from a hundred to a few thousand kms -1 . The maxi- 
mum speed, Vmax = 3750 kms -1 , is obtained for equal- 
mass maximally-spinning black holes with an t i- aligned 
spin vectors in t he orbital plane of the binary (Campan- 
el li et al.|2007b[ ). 

When Kkick is lower than the escape speed of the host, 
the MBH will wander through the galaxy in an orbit 
that depends on the detailed characteristics of the host 
potential. Several groups have calculated th e trajecto- 



ries of MBHs in smooth spherical potentials ( Madau 



Quataert 2004} |Gualandri~ 
Loeb 20081) . Asphericities i 



Merritt 2008 



2008 ) . Asphericities in the dark matter 



Blec ha & 
potential 



result m highly non-radial MBH orbits, increasing its 
decay timescale compare d to the sp herical case ( |Vicari| 
et al.|2007| |Guedes et aT7||2QQ8[ |2009[ ) . In gas-rich galaxy 
mergers, trie wandering hole may travel through a rota- 
tionally supported gaseous disk that is inhomogeneous, 
non-axisymmetric, clumpy, and turbulent. These factors 
affect the kinematics of the recoiling MBH and its abil- 
ity to accrete new fuel. In this paper, we study the or- 
bits of recoiling MBHs in high-resolution galaxy merger 
remnants that have successfully formed black hole pairs. 
These remnants are the endpoint of the 1:1 major merger 
simulation of |Mayer et al.| (| 2007[) and the 1 : 4 and 1:10 
minor merger simulations of [Uallegari et al. ([2009). We 
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construct a semi-analytical model based on the host po- 
tential formed after the merger and calculate the trajec- 
tories of MBHs for arbitrary recoil velocities and inclina- 
tion angles, accounting for the effect of dynamical friction 
and gas drag. This allows us to obtain the characteristic 
apocenter distances and return times for MBHs recoiling 
in different environments. In addition, in order to study 
the possible effect of the moving hole on the ambient gas, 
we perform high-resolution TV-body + SPH simulations 
of recoiling MBHs in the major merger galaxy remnant. 
Under the assumption that the MBH binary is able to co- 
alesce efficiently after reaching the minimum length scale 
resolved by the simulation, the hole is given a kick and 
launched into orbit in an environment characterized by 
a turbulent, actively accreting gaseous disk surrounded 
by unrelaxed stellar and dark matter components. 

The paper is organized as follows. The probability dis- 
tribution of kick velocities is reviewed in § [2] We de- 
scribe the main properties of our three host galaxies in 
§ |3j and our semi-analytical model in § [5] We assess 
the response of the gaseous disk to the motion of the 
MBH, together with its detect ability as a velocity-offset 
or spatially-offset AGN in § [6j Finally, we discuss and 
summarize our results in § [7J 

2. DISTRIBUTIONS OF KICK VELOCITIES 

The recoil velocity depends on the binary mass ratio, 
q = rai/777,2 < 1, on the magnitude of the black hole spin 

vectors = £1,2/^1,2, and on their orientation with 
respect to the orbital angular momentum vector. Fit- 
ting formulae for Vkick to fully general relativistic numer- 
ical calculations, given a particular binary co nfiguration, 



have been provided by several authors (e.g. Campanelli 



et al.||2QQ7b| |Baker etal. ] | 2QQ8| |van Me^e^t"aT||2UlU[ ). 



in the formulation of van Meter et al. ( |2U10| ): 



Vkick = ^_L,m4 + v±, s (cos£e x + sin£e 2/ ) + v\\ e z , (1) 

(2) 



V±,r. 



-H 



77^ 



V\\ = • 



(1 + <Z) 
K2V 2 + K 3 r] 3 



(0,2 cosa2 — qai cosai), (3) 
[q a\ sin a\ cos(0i — $1) (4) 



1 + q 

-a 2 sina 2 cos(0 2 - $2)] 

K s (q - l)rj 2 . 2 . , . s 

H 73— — — [q aismaicos(0i - $1) 

+a 2 sina 2 cos(</> 2 - $2)], 

where v± jm is the contribution to the recoil velocity 
from mass asymmetry, the indices || and _L refer to pro- 
jections parallel and perpendicular to the orbital an- 
gular momentum, and e x and e y are orthogonal unit 
vectors in the orbital plane. Here ai, 2 are the mag- 
nitudes of the hole spin vectors, ol\^ the angles be- 
tween ai } 2 and the orbital angular momentum vector L, 
r] = q/(l + q) 2 is the symmetric mass ratio, 0i ?2 are 
angles between the projection of along the orbital 
plane and a reference angle, and $12 are constants that 
depend on the mass ratio. Figure [l| shows a schematic 
of the binary geometry prior to coalescence. The best- 
fit parameters are A = 1.35 x 10 4 kms _1 , B = —1.48, 



H = 7,540kms-\ f = 215°, K 2 = 3.21 x 10 4 kms 
Kz = 1.09 x lQ^m s" 1 , and K s = 1.54 x 10 4 kms" 1 
Meter et al.l (12010b. 



van 




mi 



Fig. 1. — Geometry of the massive black hole binary prior to 
coalescence. 



To obtain the recoil velocity distribution as a func- 
tion of q and the probability distribution function at a 
given q, we sampl e the equation above using a method 
similar to that of |Tanaka fc Haiman] ( |2009| ) . For each 
value of the binary mass ratio we carry out 10 6 realiza- 
tions where the black hole spins have random orienta- 
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Fig. 2. — Distribution of recoil velocities as a function of binary 
mass ratio. Left panel: randomly distributed spin orientations. 
Right panel: aligned spins oriented along the orbital angular mo- 
mentum vector. In both cases the spin magnitudes are randomly 
sampled from an uniform distribution in the range < < 1. 
The solid line represents the mean of the distribution, the dashed 
lines show lcr deviations from the mean, and the dash-dotted line 
shows the maximum value attainable. The component vu vanishes 
in the case 0:1,2 = 0. 
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the distribution br oadens at higher mass ratios (see also 
Lousto et aL|2010 ). 
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Fig. 3. — Probability distribution function (PDF) of total recoil 
velocities and their components parallel and perpendicular to the 
angular momentum vector. Left panels: randomly distributed spin 
orientations. Right panels: spin aligned along the orbital angular 
momentum vector. Each panel shows the PDF for mass ratios 
q = 0.1 (solid line), q = 0.25 (dashed line), q = 0.33 (dotted line) 
and q = 1 (dash-dotted line). For 011,2 = 0, v\\ vanishes and the 
PDF is null. 

tions (with < < tt) or are aligned with the angular 
momentum vector (a i? 2 = 0). The angle — $1,2 and 
spin magnitudes are randomly sampled from a uniform 
distribution in the range — 2tt and — 1, respectively. 

Figure [2| shows the distributions of Vkick, and v\\ in 
the case where the black hole spins are either randomly 
oriented or aligned parallel to the the angular momentum 
vector. The maximum recoil velocity, ~ 3, 700 km s -1 , is 
attained in the direction parallel to L for maximally spin- 
ning black holes with anti-aligned spins oriented along 
the orbital plane, and q = 1. The mean Vkick exceeds 
100 kms -1 for mass ratios q > 0.1. If spins are randomly 
oriented, then the maximum recoil velocity attainable in 
the orbital plane is smaller than the component along 
the orbital angular momentum. The maximum kick in 
the case = is of the order of 200 kms -1 and is 
oriented in the orbital plane (y\\ = 0). 

Figure [3] shows the probability density distributions 
(PDF) 01 Kkick, v\\ , and v± at a given mass ratio, in 
cases where the spins are randomly oriented or aligned 
with the orbital angular momentum vector. The PDFs 
are calculated for mass ratios g=0.1, 0.33, 0.25, and 1. 
For randomly oriented spins, the peak of the Vkick PDF 
shifts towards lower velocities for lower mass ratios, and 



3. GALAXY MERGER REMNANTS 

We study the orbits of recoiling MBHs in three galaxy 
remnants that were previously simulated at high resolu- 
tion, and in which the MBHs of the progenitor galaxies 
were shown to successfully formed close pairs. The rem- 
nants are th e endpoint of the 1 :1 gas-rich major merger 
described in |Mayer etaT] d20Q7|), and of the 1:4 an d 1:10 
gas-rich minor mergers 01 |Callegari et al.| ( [2009 1). Fig- 
ures [4] and [5] show a 2-D rendering oi the hosts' projected 
density and the rotation curves of all components - gas, 
stars, and dark matter - of the three remnants. 

3.1. Major merger 

Our major merger remnant is p a rt of t he suit of sim- 
ulations carried out |Mayer et al.| (2007). These simu- 
lations tracked the formation ol a MBH binary down 
to parsec scales following the collision between two spi- 
ral galaxies, and were run using the massively paral- 
lel TV-body + smo othed particle hydrod ynamics (SPH) 
code GASOLINE flWadsley et 1I]|2QQ4| ). Initially, two 
equal mass galaxies composed ol 1 spherical dark mat- 
ter halo, an exponential disk, a spherical stellar bulge, 
and a central MBH, were placed on a parabolic, copla- 
nar, and prograde-prograde orbit. The galaxies have a 
virial mass of 1O 12 M and a baryonic disk mass of about 
6 x 1O 1O M , with 10% of the disk mass in a cold g as com- 
ponent (T — 2 x 10 4 if) and the rest in stars (see Mayer 
et al.|2007| for details on the structural parameters of the 
models). 

Radiative cooling and star formation were imple- 
mented until the late stages of the simulation, when the 
dark matter halos had nearly merged and the baryonic 
cores were separate d by about 6 kpc. At this point, in 
the reference run of Mayer et al. (2007), the inner 30 kpc 
of the simulation volume was refined with the technique 
of particle splitting to achieve a spatial resolution of 2 
pc. The star and dark matter particles were not splitted 
to limit the computational burden. The calculation that 
we use in this paper has a softening of 10 pc and a mass 
resolution in the gas component of 5000 M after split- 
ting (for a total of N g = 1.4 x 10 6 SPH particles). We 
use this run rather than the hig hest resolution simulation 
presented in Mayer et al. (2007|) because this resolution 



matches those ol the minor merger runs also used in this 
work, and because it allows to perform more recoil exper- 
iment s at a reasonable computational cost ( Mayer et al. 
2007 verified the convergence in the mass distribution ol 
the merger remnants between the different resolutions). 
The radiation physics in the refined region was modeled 
via an effective equation of state, an ideal gas with adia- 
batic index 7 = 7/5. Calculations that include radiative 
transfer have shown that this simple treatment approx- 
imates well the thermodynamics of a solar metallicity 
gas heated by a star burst over a wide range of densities 
flSpaans fc s5k|20QQ| |Klessen et aT|2007[ ). 

By the end ol the simulation, t = 5 Gyr after the 
merger, tidal tails continue to evolve in the outer parts 
of the remnant and over 60% of the gas has been fun- 
neled to the central region, forming a compact nuclear 
gaseous disk of mass = 3 x 1O 9 M , radial length 
a = 75 pc, and scale height b = 20 pc. The inner disk 
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is turbulent and rotationally supported. The two MBHs 
have sunk down from about 40 pc to a few parsecs in less 
than a million years, are gravitationally bound to each 
other, and their orbital decay is controlled by gas drag in 
the supersonic regi me, not dynamical friction against the 
stellar backgro und flMayer et al.|2007| |Escala et al.|2004| 
Ostriker 1999). Figure [5] (left panel) shows the rotation 
curve tor the different components of the galaxy remnant 
at the endpoint of the merger. The gaseous component 
dominates the dynamics of the inner 300 pc of the re- 
mant galaxy. The escape speeds from the center in the 
perpendicular (9 = 90°) and parallel (9 = 0°) directions 
to the angular momentum vector of the gaseous disk are 

= 1,400 kms -1 and v\ = 1,320 kms -1 , respectively. 

3.2. Minor mergers 
The hierarchical ACDM structure formation scenario 



predicts a high frequency of minor gala x y mergers ( Lacey 



fc Cole 1993 Fakhouri fc Ma 2008). Callegari "etaT 



( |2UU9p run a set of collisionless^Mry" , with f g = D) 
and gasdynamical ("wet") minor mergers with the same 
gas fraction in the primary and secondary galaxies, either 
f g = 0.1 or 0.3. As in the major mergers, the primary 
galaxy is a Milky Way-mass galaxy sampled with 10 6 
particles dark matter particles and (initially) 2 x 10 5 star 
and 10 5 gas particles. The secondary galaxy is obtained 
by rescaling the primary galaxy in mass by 1/ 4 or 1/10 
of its mass, according to the scaling laws of |Mo et al. 
( |1998| ) for CDM galaxies, and is sampled by an equiv- 
alent number of particles as the primary. In the 1:10 
merger, which is assumed to take place at z = 3 since 
this is an ideal redshift window for LISA, the masses, 
radii and velocities of the galaxies are s caled further by 
a fact or H(z)/Hq as expected in CDM ( |Callegari et al. 
20091 ). 

While no MBH pair was formed in the collisionless runs 
(the two holes were defined as a "pair" if their relative 
orbit shrunk down to a separation equal to twice the soft- 
ening, which sets the resolution limit), the presence of a 
gaseous component allows the formation of a MBH pair 
because the massive gaseous cores in which the MBHs 
are embedded sink efficiently by dynamical friction de- 
livering the holes at the center of the remnant. Further- 
more, the higher background density arising as inflows 
driven by tidal torques associated with bar instabilities 
and shocks push gas to the center of the galaxies during 
the merger and enhance dyn amical friction in t he late 
stages of the orbital decay (Calleg ari et al.||2009 ). Here 
we consider two such dissipational runs: a 1:4 merger 
with initial gas fraction f g = 0.1, and a 1:10 merger with 
f g = 0.3. Both resulted in the formation of MBH pairs 
separated by twice the softening length (200 pc for the 
1:4 merger and 40 pc for 1:10), and it is expected that 
the holes will subsequently m erge on a timescale < 1 
Gyr (see |Callegari et a l. 2009). The galaxies used in the 
simulation were initialized with a mixture of primordial 
metallicities. The runs were performed with the GASO- 
LINE code, including this time radiative cooling, atomic 
line cooling (by H and He), and cooling via collisions, star 
formation, and blast-wave supernova feedback until the 
end of the simulation, rather than switching to an effec- 
tive equation of state. The resulting interstellar medium 
(ISM) density distribution is clumpier than in the case 




Fig. 4. — Projected gas density of the galaxies hosting a recoiling 
MBH. (a) Remn ant of a 1:1 merg er simulation with initial gas 
fraction f g = 0.1 ( |Mayer et al.|2007| . (b) Re mnant of a 1:4 merge r 
simulation with initial gas traction j g = 0.1 JCallegari et al. 12009). 
(c) Remn ant of a 1:10 merger simulation with initial gas traction 
f g = 0.3 (Call egari et al.||2009[ ). The scaling is logarithmic and 
the highest and lowest values of the scale are given in units of 
atoms cm -3 . Note that the lower central density in the 1:4 merger 
remnant, which manifests as a shallower profile of the inner rotation 
curve inside a few hundred parsecs relative to the 1:10 merger case, 
is an artifact of the larger gravitational softening adopted in this 
run. 

of the major mergers that adopt an effective equation of 
state and neglect supernova explosions. 

The most important difference between the major and 
the minor mergers, however, is that bar-driven torques 
are much weaker in the latter case. This is especially 
true in the primary galaxy, which is only weakly destabi- 
lized by the secondary: as a result, gas inflows are much 
weaker, and < 10% of the gas that was originally at 
large radii reaches the central hundred parsec region by 
the end of the merger, as opposed to more than 60% in 
the case of major mergers. Figure [5] (middle and right 
panels) shows the rotation curves of the minor merger 
remnants. They include gas that will eventually form 
stars; furthermore, heating and stirring by supernova 
feedback will make the gas more pressure-supported and 
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Fig. 5. — Rotation curves of the galaxy remnants hosting a recoiling MBH. Left 'panel: 1:1 merger. Middle panel: 1:4 merger. Right panel: 
1:10. The colors represent the circular velocity profiles of the gas (red), stars (blue), dark matter (green), and total (black), respectively. 



extended. Therefore, we expect the MBH return time 
in minor mergers to be largely determined by dynamical 
friction against the stars and the dark matter. In ad- 
dition, the lower central densities imply a lower escape 
speed, so we expect that holes with a given kick veloc- 
ity will travel larger distances in this case compared to a 
major merger. The escape speeds from the center in the 
perpendicular [6 = 90°) and parallel (0 = 0°) directions 
to the angular momentum vector of the gaseous disk are 

= 520 km s -1 and v\ = 490 km s -1 in the 1:4 merger, 



respectively, and vjr 
in the 1:10 merger. 



490 kms and vl = 450 km s 



4. SIMULATIONS OF RECOILING MBHS 

We performed six si mulations of recoiling MBHs in the 
1:1 merger remnant of |Mayer et al.| ( |2007| ). A MBH bi- 
nary of total mass M. = 5.2 x 10 b M was replaced with 
a single black hole of the same mass. The hole was given 
a kick of Vkick = 400, 800, and 1,200 kms -1 at an an- 
gle = 0° or 90° relative to the angular momentum 
of the disk. In Section [2] we showed that maximum re- 
coil velocities are of order 200 kms -1 if spins are aligned 
with the orbital angular momentum of the binary. Here, 
we make no assumptions regarding the spin orientation 
or the orientation of the binary orbital plane relative to 
the gaseous disk, since the main objective of the simula- 
tions is to calibrate the semi-analytical model described 
in Section [5] However, these configurations will enter in 
the calculation of the probability of observing recoiling 
MBHs as off-nuclear AGN during their wandering phase 
(see Section [6J. For comparison, we have also run a 
simulation without the MBH. 

Our gas mass resolution is 5 x 1O 3 M , for a total of 
1.5 x 10 6 SPH particles, and the gravitational softening 
is 10 pc for the MBH as well as the SPH particles. The 
hole was kicked starting from where the binary was, at 
the center of mass of the gas distribution located some 
20 pc away from the center of the potential. The simu- 
lations were run using GASOLINE, and were restarted 
using the same the rmodynamica l para meters as in the 
refined part of the Mayer et al. ( 2007[ ) simulation, i.e. 
the SPH particles were treated as an ideal gas with adi- 
abatic index 7 = 7/5. Since only a small fraction of the 
nuclear gas turns into stars during the timescale of the 



simulations, we neglect star formation to reduce the com- 
putational burden. We run the simulations for up to 50 
Myr, a compromise between obtaining a good calibration 
for our semi- analyt ical models and lim it CPU time. It 
has been shown by |Dotti et ah] ( |2007[ ) that MBHs mov- 
ing in nuclear disks with varying mass fraction of stars 
produce a drag (via dynamical friction) that is almost 
indistinguishable from the purely gaseous case. It is gas 
thermodynamics instead that determines the magnitude 
of frictional drag in both purely gaseous or star forming 
disks, as it sets the density and temperature distribution 
of the gaseous disk ( |Escala et al.|2005| |Mayer et al. 2007 ) . 



It is important to emphasize that, in a model with an 
effective equation of state, the temperature should not 
be associated with the physical temperature of a spe- 
cific phase of the interstellar medium in the nuclear disk, 
e.g. t he molecular gas that m akes up most of the disk 
mass (Downes & Solomon 1998), but rather with an effec- 
tive pressure that provides a realistic match to the pres- 
sure scale height of the disk. In a realistic case the cold 
molecular phase is embedded in a hot ionized medium 
that fills most of the volume and sets the actual pressure 
scale height. The pressure scale height of our nuclear 
disk is zo =15 pc at the beginning of the simulations, 
and Zq =21 pc at t =20 Myr, and compares well with 
the pressure scale height of 23 pc observed i n Mrk 231 
( |Downes fc Solomon) 1998| |Davies et al.||2004| ). Likewise, 
the effective sound speed used throughout the paper (cal- 
culated directly from the simulations) should be close to 
the sound speed of the hot, diffuse phase (the cold molec- 
ular phase has a sound speed of < lkms -1 ). The peak 
of the sound speed in the major merger is 200 kms -1 , 
while it is 34 kms -1 and 14 kms -1 in the 1:4 and 1:10 
mergers, respectively. 

Figure [6] shows the resulting orbits for MBHs kicked 
with Vkick = 800 kms -1 and Vkick = 1,200 kms -1 in a 
direction that is either along the disk plane (cases a and 
c) or perpendicular to it (cases b and d). The boxes have 
varying lengths and the orbits are centered on the point 
where the potential is the deepest. The orbit for the 
case Vkick = 1, 200 kms -1 and = 90° is shown again in 
Figure^ together with the evolution of the background 
gas number density and specific internal energy. 
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Fig. 6. Orbit of a MBH of mass M. = 5.2 x 10 6 M in major 
merger remnant. The panels show a hole kicked with (a) Vkick — 
800 km s -1 and = 90° (parallel to the gaseous disk plane), (b) 
Hick = 800 kms" 1 and = 0°, (c) V kick = 1, 200 kms" 1 and 9 = 
90°, and (d) V kick = 1, 200 kms" 1 and 6 = 0°. The colors indicate 
time along the orbit, from t < t r un/3 (brown) to t > 2£ run /3 
(yellow). The runtime of the simulation is t run = 10 Myr for 
Hick = 800 kms -1 and t run = 20 Myr for V kick = 1,200 kms -1 . 
The black dots mark 1 Myr intervals. Orbits are plotted with 
respect to the deepest point in the disk potential. 

5. SEMI-ANALYTICAL MODEL 

To obtain a large sample of MBH orbits at relatively 
low computational cost, we have constructed a semi- 
analytical model based on our numerical results. We can 
then follow the trajectories of MBHs until they return to 



the center of their host, on timescales that can be much 
longer than the time it is practical to run a full TV-body 
+ SPH simulation, particularly for larger kick velocities. 
The motion of a MBH in a multi-component potential 
can be studied analytically by integrating numerically 
its equation of motion under the action of a conservative 
force V<3> and a damping frictional term: 



dv 
~di 



-V$ + /df, 



(5) 



where the total potential is $ = <I>g as + $ stars + $dark- 
The gradient of the potential is computed directly from 
the simulation data. The dynamical friction term for 
an object of mass M. moving with relative velocity v 
through a gaseous medium of density p is given by 



DF 



47rG 2 M.lnAp 



(6) 



The velocity integral, I v = I(M), depends on the Mach 
number M = v/c s , where c s = ^^ksT/p is the sound 
speed. In our model, the local gas density and sound 
speed are measured directly from the simulation data in 
a sphere of radius 



#sphO) = # sp h(lkpc) 



p(lkpc) 



1/3 



(7) 



around the MBH, where p(r) is the average density at 
a distance r from the host center and we choose a ref- 
erence radius i? sp h(lkpc) = 100 pc. As for the shape of 
the velocity integra l, we a dopt the hybrid formulation of 



Tanaka fc Haiman l (|2009[ ), who combined the p rescrip- 
tions of [Ostriker |T9gUi rand |Escala et aT] ( |2004[ ) into: 



I(M) 



- J lM exp 



(-*): 



0.5 In A erf(^| 
(0 < M < 0.8); ' 
1.5 In A [erf (^) - ^Mexp (-*f 
(0.8 < M < 1.5); 
0.5 In (1 - jfo) +lnA 
{(M > 1.5). 



(8) 



The C oulomb logarithm is set to In A = 3.1 |Escala et al. 
'2004). In a collisionless medium, the velocity integral 



epends on the local velocity distribution of the stars or 
dark matter. If the medium is isotropic and the velocity 
distribution is Maxwellian, then I v = I(X = v/y/2a), 
where a is the one-dimensional isotropic velocity disper- 
sion. Equation ([6| then b ecomes the Chandras ekhar dy- 
namical friction formula ( Ch andrasekhar|[l943| ) , where 



I(X) = erf (X) 



2X 



(9) 



and p is the stellar or dark matter density. In a recent 
merger, however, the velocity distribution of the stellar 
and dark component will deviate from Maxwellian. We 
therefore calculate the velocity integral numerically from 
the actual velocity space distribution of the particles in 
the simulations, 



JO 



f(v')dv' 



(10) 



7 




Fig. 7. The orbit of a recoiling MBH with M m = 5.2 x 10 6 M©, y kick = 1,200 kms" 1 , and = 90° (black dots), as in panel (c) of 
Figure 6. The white arrows show the instantaneous direction of the velocity vector at the time of the corresponding snapshot. The color 
coding indicates gas number density in units of log(n/cm — 3 ) (left panel) and specific internal energy in units of log(it/ km 2 s — 2 ) (right 
panel). 

where f(v')dv' is the number of particles with velocities using an adaptive Adams-Bashforth-Moulton algorithm, 
between v' and v' + dv' in the vicinity of the hole. The Figure [8] shows a comparison between the orbits of 

equation of motion of the MBH is integrated numerically a MBH in our numerical simulations versus the semi- 
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Fig. 8. — Comparison between the orbit of a recoling MBH of mass M m = 5.2 x 10 6 Mq in our numerical simulation versus the semi- 
analytical model, (a) Same simulation as in panel (a) of Figure 3. (b) Semi-analytical orbit. For clarity, panels (a) and (b) show the orbit 
10 Myr after the initial kick, (c) Radial evolution of the orbit in the semi-analytical model (gray) compared to the simulation (red). 



analytical model. The hole was kicked with Vkick = 
800 krns" 1 in the gaseous disk of the major merger rem- 
nant. The trajectory is shown in projection in the left 
and middle panels, while the right panel depicts the ra- 
dial evolution of the orbit. Over the timescale covered 
by the simulation, the semi-analytical model is clearly 
able to reproduce both the amplitude of the oscillations 
(dictated by the potential) and the decay by friction. 
Figure [9] shows the first apocenter distance R max ^ esti- 
mated from our semi- analytical approach, for MBHs re- 
coiling in our 1:1, 1:4, and 1:10 merger remnants along 
a direction that is either parallel or perpendicular to the 
disk angular momentum vector. The black hole mass 
was scaled in each case to match the M % — <j relation, 
yielding M. = {5.2,3.2,2.8} x 10 6 M for the 1:1, 1:4, 
1:10 mergers, respectively. At supersonic speeds, gas 
drag plays an important role in damping the MBH or- 
bit in the 1:1 remnant. Shortly after the kick, when 
the hole's is moving at Mach number 1 < A4 < 5, gas 
drag accounts for over 70% of the total frictional force. 
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Fig. 9. — First apocenter distance of recoiling MBHs as a function 
of kick velocity. Red: 1:1 merger remnant. Orange: 1:4. Green: 
1:10. Solid curves: MBHs recoiling in a direction perpendicular to 
the disk angular momentum vector. Dashed curves: MBHs recoil- 
ing in a direction parallel to the disk angular momentum vector. 



In the deep potential well of the remnant, the gas den- 
sity can reach 10 5 particles cm -3 , and MBHs kicked in 
the plane of the disk with initial velocities as high as 
1,200 kms -1 reach apocenter distances of only 0.3 kpc. 
As Vkick approaches the escape speed along this direction, 
i.e. Vkick = 1,350 kms -1 = 0.98v^, the first apocenter 
distance reaches i? max = 10 kpc. 

MBHs recoiling in the direction parallel to the disk 
angular momentum vector can reach the low-density 
outer regions for slightly lower initial kicks, and there- 
fore reach larger apocenter distances. A hole with Vkick = 

1,200 kms -1 = 0.91i>J reaches i? max ~ 1 kpc and has a 
return timescale that is close to 10 Gyr. Minor merg- 
ers have significantly lower central density profiles, and 
apocenter distances can be significantly larger even for 
modest recoil speeds. In the 1 : 4 merger remnant, a 
kick of Vkick = 500 kms -1 perpendicular to the disk can 
displace the hole by i? m ax ~ 10 kpc. The same initial 
speed yields a maximal separation of i? max = 100 kpc in 
the 1 : 10 case. 

6. DETECTABILITY OF RECOILING MBHS 

Three main methods have been proposed to detect 
a recoiling hole: 1) through the imprint of the moving 
MBH on the gaseous medium of the host galaxy; 2) as a 
kinematically-offset AGN, as the broad line region (BLR) 
carried by the hole is displaced by a velocity Av from the 
narrow starlight associated with the host galaxy; 3) as 
a spatially-offset AGN, where the accreting MBH is ob- 
served to be displaced by a distance Ar from the host's 
nucleus. 



6.1. Interaction with the gaseous disk 
Devecchi et al. (2009[) recently carried out a suite of 



numerical simulations aimed at studying the detectabil- 
ity of recoiling MBHs in the hot gas of elliptical galax- 
ies. The host galaxy was mod eled as a spherical Hern- 



quist sphere (Hernquist 1990), with a mass resolution 
of 2 x 1O 7 M and a polytropic equation of state with 
7 = 5/3. They found that, at high Mach numbers, the 
shocked gas can produce an X-ray signature that may be 
detectable with the Chandra satellite in nearby galax- 
ies. Simulations o f recoiling holes s urrounded by a thin 
circumbinary disk ( |Rossi et al.|2010 ) show that the likeli- 
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hood of observable emission is low and depends strongly 
on the orientation of the kick relative to the thin disk, on 
the magnitude of the kick, and on the mass of the disk. 

In the numerical simulations described in Section[4j the 
inner disk is rotationally supported (v rot = 300 km s -1 ), 
turbulent (^turb = 500kms _1 at r = 20 pc), and ac- 
tively accreting material from the outer regions. Dense 
columns of cooler gas reach the nuclear region, constantly 
changing its internal energy profile. We find that the ra- 
tio between the integrated energy deposition of the MBH 
(associated with energy losses from gas frictional drag) 
and the integrated change in internal energy by the gas 
(measured in the simulation without MBH) is very small, 
AE-pp / AE[ nt = 0.02. As a result, overdensities produced 
in the gas by the recoiling M BH are rapidly weak ened. 
In contrast with the results of |Devecchi et al ] p009), the 
removal of the central hole from the disk potential does 
not have a significant effect on the central density pro- 
file, since the mass of the hole amounts to only 2% of the 
mass of the nuclear disk. At the end of the major merger 
+ MBH simulation, the gas density and internal energy 
profiles (Fig. [7|) are nearly undistinguishable from those 
in the no-MBHrun. 

By contrast, in the 1:10 merger, the integrated change 
in gas internal energy is comparable to the orbital energy 
lost by the MBH over a few tens of orbits. High Mach 
numbers can be reached at lower recoil speeds since the 
sound speed profile rises slowly to reach 30 kms -1 at 
50 pc and then stays constant out to 100 pc (whereas 
in the major merger the sound speed reaches a peak of 
230 kms -1 at 10 pc and then drops to 150 kms -1 at 100 
pc). This suggests that minor merger remnants may be 
suitable sites for the detection of recoiling MBHs through 
their interaction with ambient material. 

6.2. Off-nuclear AGNs 

The detectability of MBHs as offset AGNs depends on 
the hole's phase-space location at the time of observa- 
tion and on its ability to accrete material and shine at 
that location. We assume that, upon election, the hole 



Loeb|[2007| ) and 



carries a punctured accretion disk (e.g. 
an associated BLR that is visible for a time comparable 
to the visc ous lifetime of the disk. The mass car ried by 



the disk is flLoeb|2007| |Volonteri fc Madau||2008l ) 

M disk = 1.9 x 10 6 a:?- 8 C- - 6 M 7 2 - 2 K- 2 - 8 M„ 



(11) 



/ 3 ivl , 

where a_i is the viscosity parameter in units of the 
fiducial value a — 0.1, Mj is the MBH mass in units 
of 1O 7 M , ( = e_i//#, e_i is the matter-to-radiation 
conversion efficiency in units of 10%, Je = M/Me is 
the hole's accretion rate in units of the Eddington rate 
M E = 4:7rGM m m p /(eaTc) (we use Je — 0.1), and Vs is 
the kick velocity in units of 10 3 kms -1 . The condition 
^disk < M. requires 



F kick > 550a:f 8 C 



-0.21 7i/f0.43 



M^ 6 krns" 



(12) 



When the above equation is not satisfied, we impose 
Mdisk = M 9 : this yields an AGN lifetime equal to 
t AGN M./M = eca T /(^Gm p f E ) ~ 45CMyr. The 
condition M^isk < M % also sets the maximum outer ra- 
dius of the disk 



i?out<150a:;- 5 V°- 43 M 7 -°- 86 i? s 



(13) 



where R s = 2GM,/c 2 is the Schwarzschild radius of the 
hole. At large radii (> 0.1 pc), gas in the disk may 
be gravitationally un stable and prone to star formation 
dCollin fc ZaTm|[l999l ). However, the maximum disk size 
considered in our model is ~ 10 -3 pc, and therefore we 
neglect this effect in our analysis. We also find that ram 
pressure stripping will not affect the accretion disk as 
the MBH travels through the host. For a hole moving 
with velocity v relative to a medium of density p, we can 
calculate the stripping radius R str by equating the ram 
pressure force F ram = pv 2 to t he restoring force of th e 
disk, using the density profile of Goodman & Tan (2004). 
In the cases considered here, i?o U t <C i4tr, and therefore 
the disk is able to survive. 

Because precise measurements of the peaks of broad 
emission lines in quasar spectra are challenging, the ob- 
servation of potential recoil candidates as kinemtically- 
offset AGNs - whose broad emission lines are substan- 
tially shifted in velocity relative to the narrow-line gas 
left behind ([Bonning et al.|2007 ) - is biased toward large 
kicks. In the following we shall assume that a recoil- 
ing hole is detectable as a kinematically-offset AGN if 
its velocity shift, Av, is above the detection threshold 
Avt = 600 kms -1 . At every time step of our numerical 
integration of the hole's orbit, we calculate the instan- 
taneous accretion rate, the associated luminosity, and 
update the MBH mass. When the initial bound disk is 
exhausted, the BLR is assumed to vanish and the recoil- 
ing MBH is no longer recognizable as such from kine- 
matics measurements. From this point on, the MBH is 
allowed to accrete material from the ISM. The fueling 
rate depends on the local density and sound speed of the 
gaseous medium and on the speed of the hole following 
the Bondi-Hoyle formula: 



M B = 7 x 10" 7 n o M^ 3 " 3 M yr" 1 , 



(14) 



where no is the local gas number density of the ISM in 
units of cm -3 , vs = (v 2 + c 2 ) 1 / 2 is expressed in units 
of 10 3 kms -1 , and v is the instantaneous MBH velocity 
relative to the host. The mass accretion rate in this phase 
is given by the minimum between the Bondi rate and the 
Eddington rate, and the bolometric luminosity is then 



L = ec 2 x mm(M B ,M E ). 



(15) 



To be detected as a spatially-offset AGN, a recoiling 
accreting hole must be both above a threshold spatial 
separation Ar t from the host center as well as above a 
luminosity threshold L t . We set Ar t equal to roughly 
twice the resolving power of the telescope - 0.2" for 
the Hubble Space Telescope (HST) and the James Webb 
Space Telescope (JWST), and 1" for Chandra. For ref- 
erence, at redshift z = 1 (0.1) an offset of 10 kpc cor- 
responds to an angular separation of 1.25" (5.4"). In 
the following, we set the threshold X-ray luminosity to 
Lt.x fy^^D 2 



fii m 4iTD 2 L where f\ 

im 



3.8 x 10" 16 ergcm- 2 s -1 



is the on the on- axis Chandra limit and Dl is the lu- 
minosity distance, and use a bolometric to X-ray correc- 
tion of Y2- mkeV = 20 if L/Le < 0.1 a nd X2-iokeV = 50 
otherwise ( Vasudevan & Fabian |2007[ ) . The optical B- 
band luminosity threshold L t ^ is calculated using a lim- 
iting flux of /ii m = 2 x 10~ 14 , erg cm" 2 s -1 from observed 
broad-line AGN (Zakamska et al.||2003|) and the B-band 



bolometric correction is taken from Marconi et al. (2004) 
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Both X-ray and optical thresholds are consistent with ob- 
serv ational limits from mult i- wavelength studies of AGN 
(e.g pierce et al.||2010 ). 



1200 



to 

E 



800 



400 



Bondi-Hoyle- 

Lyttleton 

accretion 




Accretion 
from initial 
disk, BLR 



Fig. 10. — Phase diagram of the orbit of a MBH of mass 
M m = 5.2 x 1O 6 M recoiling with V kick = 1,300 kms -1 along 
the nuclear disk of the major merger remnant at z = 0.025. The 
MBH undergoes disk accretion and is shining with an X-ray lu- 
minosity > L t x along the thin red line, and is powered by ISM 
accretion in the orange region. 



Figure [To] shows a breakdown of the orbital phase 
diagram ofa MBH of mass M. = 5.2 x 10 6 M recoiling 
with velocity Vkick = 1,300 kms -1 in the plane of the 
disk of the major merger remnant at z = 0.025. Fiducial 
regions where Av > 600 kms -1 , Av > 1,000 kms -1 , 
and Ar > 2" are marked in the velocity-distance plane. 
In this example 2" correspond to 1 kpc. The MBH 
undergoes disk accretion and is shining with an X-ray 
luminosity above L t x along the thin red line: it is 
therefore detectable as a kinematically offset AGN when 
it is still close to the center, along the portion of the 
red line where Av > 600 kms -1 . The MBH is instead 
detectable as a spatially-offset AGN in the low-velocity 
region where it accretes rapidly from the ISM (orange) 
or from the initial disk (red) and where Ar > 1 kpc at 
this redshift. Note that large kinematic offsets occur 
near pericenter and can be observable for a shorter 
timescale than large spatial offsets, which occur near 
apocenter. In the following, we study in details the 
detect ability of recoiling MBHs in two cases. 

Case 1 . A best-case scenario for the detectability of 
recoiling MBHs as offset AGNs, in which: (a) the 
spins are randomly oriented with respect to the orbital 
angular momentum prior to the merger of the MBH 
binary (a possible mechanis m for misalignment o f the 
spins is chaotic accretion, see King & Pringle 2006 k; (b) 
the orbital plane of the binary is randomly oriented with 
respect to the plane of the gaseous disk (e.g. because 
of turbulent motions at sub-pc scales that may be 
uncorrelated with the angular momentum of the gaseous 



disk); and (c) 1:1, 1:4, and 1:10 galaxy merger remnants 
host MBH binaries of mass ra tios q = 1,0.2 5 ,0.33 , 
respectively. In a recent study, |Callegari et al. (2010) 
showed that during a 1:10 merger the primary galaxy 
may induce strong tidal torques onto the satellite that 
can result in the rapid growth of the smaller black hole, 
yielding in some cases binary mass ratios of q = 0.5. In 
simulations where the gas fraction was set to f g — 0.3 
and the merger was coplanar, a setup similar to the our 
initial 1:10 merger, the final mass ratio of the binary 
was found to be q = 0.33. 



Case 2 . A worse-case scenario for the detectability of 
recoiling MBHs as offset AGN, in which: (a) the spins 
are aligned with the orbital angular momentum vector of 
the binary prior to the merger event. Studies of the spin 
evolution of MBHs in a gaseous medium show that spins 
tend to align in gaseous environments due to frame drag- 
ging, yielding maximum re coil velocities of ~ 200 kms -1 
(Bogdanovic et al. 2007). This effect drastically de- 
creases the probability of observing both kinematically 
and spatially offset AGNs in gaseous environments; (b) 
the orbital plane of the binary is aligned with the gaseous 
disk. This assumption, together with (a), implies that all 
recoils occur in the plane of the disk; and (c) 1:1, 1:4, and 
1:10 galaxy merger remnants host MBH binaries of mass 
ratios q = 1, 0.25, 0.1. 

In both cases we consider a random uniform distri- 
bution of spin magnitudes in the range < < 1. 
Conservation of linear momentum also requires that the 
mass of the disk carried by the hole reduces the initial 
kick velocity to an effective value 



VeS = Vkick 



M. 



M. + M di 



disk 



(16) 



We account for this effect when estimating the detection 
probabilities associated with recoiling MBHs below. 




200 400 600 800 1000 1200 1400 200 400 600 800 1000 1200 1400 
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Fig. 11. — Relevant timescales for assessing the detectability of 
recoiling MBHs as kinematically-offset AGNs in the major merger 
remnant, as a function of the kick velocity. Gray curves: re- 
turn timescale of a MBH of mass M. = 5.2 x 1O 6 M . Dashed 
cyan curves: total time during which the MBH has velocity offset 
Av > 600 kms -1 , i.e. above the threshold for the detectability 
of kinematically offset AGNs. Solid cyan curves: total time over 
which the MBH has a velocity offset Av > 600 kms -1 while shin- 
ing with an optical luminosity > L tj s (for z = 0.1). The dashed 
blue and solid blue curves use a threshold for detection of kine- 
matically offset AGNs of Av > 1,000 kms -1 . Solid black curve: 
lifetime of the accretion disk carried by the hole as a function of 
recoil velocity. 
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6.3. Kinematically- offset AGNs 

Several candidates of velocit y-offset AGNs have been 
discussed in the literature (e.g. i Magain et al. 2005; [Ko-1 
mossa et al.|2008||Shields et al.|2UU9[ |Civano et al.|2U10p , 



but none has b een positively c onfirmed to date. 

The study of |Bonning et al. (2007) shows that there is 
only scarce (if any) evidence tor recoiling MBHs in Sloan 
Digital Sky Survey quasars. 

To assess the detectability of recoiling MBHs as kine- 
matically offset AGNs, we have generated a library of 
orbits using the semi-analytical model described in Sec- 
tion [5j for kick velocities lower than the escape speed 
from the host and return timescales shorter than the 
Hubble time. We have computed the accretion rate and 
the resulting optical and X-ray luminosities along the 
trajectory of the recoiling holes. 

Figure [TT] shows the return timescales of recoiling 
MBHs in tne major merger remnant, for different kick ve- 
locities in the plane of the nuclear disk and perpendicular 
to it. The figure also shows the lifetime of the initial disk 
carried by the hole, the total time spent by the hole above 
the threshold velocity offsets Av t > 600 kms -1 and 
Avt > 1, 000 kms -1 , as well as the detectability window 
as a velocity-offset AGN, i.e. the total time during which 
the MBH is both above the velocity and the luminosity 
thresholds (at z = 0.1). The detectability window of a re- 
coiling hole depends on the kinematic profile of its orbit. 
For moderate kick velocities, where the MBH can un- 
dergo many pericenter passages before the initial disk is 
exhausted, more eccentric orbits yield higher relative ve- 
locities, making the MBH observable as a kinematically- 
offset AGN. When the MBH is kicked in the direction 
perpendicular to the disk, the return timescale is longer 
and orbits are less eccentric, yielding generally shorter 
observability windows for large kinematic displacements. 
In our major merger, a MBH launched with recoil veloc- 
ity Vki ck = 1,200 kms -1 in the = 90° direction, will 
return to the center of the galaxy in about 0.3 Gyr, and 
would be observable as a kinematically-offset AGN for 
only 0.2 Myr. The detectability window decreases by a 
factor of 4 in the = 0° case for the same recoil speed. 
In the minor mergers the threshold velocity is reached 
only for kicks with Vkick > v esc , which are are not shown 
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Fig. 12. — Probability of observing a recoiling MBH in the ma- 
jor merger remnant as a kinematically-offset AGN. Cyan curves: 
probability that the MBH has a velocity offset along the line of 
sight A^i os > 600 kms -1 while shining with an optical luminosity 
> L t ^B (for z = 0.1). Blue curves: same for Ai;i os > 1, 000 kms -1 . 
Both probabilities were computed assuming our best-case scenario 
Case 1 (see text for details). 



here. 

The probability of observing a recoiling hole as a 
velocity-offset AGN with Av\ os along the line of sight 
(LOS) can be written as 

^Avi os (Vkick) = P{Avi os >Av t ,L>L t } (V kick ,z) (17) 
xP V kick (ai,a2,^), 

where P{Av los >Av t ,L>L t }(V kick ) is the probability that 
a hole with initial kick Vkick has an velocity offset 
A^ios > Av t while moving along the LOS at the time 
of observation and shining with bolometric luminosity 
above the threshold L t (assuming no obscuration). Here, 
Py kick (ai, c?2, q) is the probability that a coalescing MBH 
binary with spin vectors d[ , a~2 and mass ratio q recoils 
with velocity Vkick > Av t . In terms of the relevant 
timescales, 



PAv los (Vkick, z) = 



t{Av los >v t ,L>L t } 

(Vk ick ) 



^return 
X P V kick («l,«2,g), 



(18) 



where t {Avios > VuL > Lt y(V kick , z) is the integrated time 
during which A^i os > Av t alon g th e LOS, and L > L t . 

We have evaluated equation (18) for a range of recoil 
velocities, converting each value oi Vkick into an effective 
recoil velocity V G s (using Equation [l6| ), and interpolat- 
ing the values of PfAvi os >Av t ,L>L t } using the timescales 
shown in Figure [ll] The resulting probability is shown 
in Figure [12] under the assumptions of our best-case sce- 
nario Case 1 (random alignment of the binary orbital 
plane and spin). The probability peaks at ~ 800 kms -1 
to a maximum value of only Pavi os ~ 5 x 10 - 7 for 
A^ios > 600 kms -1 , and is nearly 2 orders of magni- 
tude smaller for A^i os > 1,000 kms -1 . We calculated 
the same probability in Case 2, where the spins of the 
binary are aligned with the orbital angular momentum 
vector, the orbital plane is aligned with the gaseous disk, 
and all kicks occur in the plane of the disk. In this sce- 
nario the maximum recoil velocity is Vkick ~ 200 kms -1 
and the holes are unobservable as kinematically-offset 
AGNs. 

6.4. Spatially- offset AGNs 

We can use a similar method to that described above to 
estimate the observability of recoling MBHs as spatially- 
offest AGNs. Figure 13] shows the integrated time during 
which the MBH is displaced by a threshold distance of 
Ar t = 0.4" and Ar t = 2" at z = 0.1 (corresponding 
to twice the resolution of HST/JWST and Chandra, re- 
spectively), as well as the integrated time during which 
the MBH is both displaced by Ar t and is shining with 
L > L t . The upper, middle, and lower panels show re- 
sults for the 1:1, 1:4, and 1:10 mergers, respectively. 

In the major merger, only large recoil velocities (of the 
order the escape speed) can displace the MBH by Art 
when the kick occurs in the plane of the gaseous disk. 
The hole can shine for timescales far longer than the 
lifetime of the accretion disk carried along, since there is 
a large supply of gas in the galactic disk for fueling near 
apocenter, where the MBH spends most of its wandering 
time. 
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Fig. 13. — Relevant timescales for assessing the detectability of recoiling MBHs as spatially-offset AGNs at z = 0.1, as a function of the 
kick velocity. Gray curves: return timescale of a MBH recoiling in the 1:1, 1:14, and 1:10 merger remnants. The displacement threshold is 
set to Art = 0.4". Dashed green curves: total time during which the MBH has a spatial offsets of Ar > Art = 0.4". Solid green curves: 
total time during which the MBH has a spatial offsets of Ar > Art = 0.4" while shinning with luminosity L > L t ,B- Dashed brown lines 
and solid brown lines: same total times but for the thresholds Art = 2" and L > L t x- Solid black line: the lifetime of the initial disk as 
a function of recoil velocity. 
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Fig. 14. — Probability of observing a recoiling MBH as a 
spatially-offset AGN. The colors represent the 1:1 (blue), 1:4 (or- 
ange), and 1:10 (red) merger remnants. Solid lines: detectability 
thresholds of Art = 0.4" and L > L t ^B- Dashed lines: Art = 2" 
and L > L tj x- All probabilities were computed assuming our best- 
case scenario Case 1 (see text for details). 



In the minor mergers, the potential wells are shallower 
and large apocenter distances can be reached at low kick 
velocities (Vkick < 500 kms -1 ). 

For trajectories having low inclinations relative to the 
gaseous disk, Bondi accretion can be significant and the 
MBH can wander sufficiently far from the center of the 
host to be detectable as an offset X-ray AGN at z < 0.1. 
The detectability window typically ranges from a few to 
tens of Myr. 

The probability of observing a recoiling MBH as a 
spatially-offset AGN with displacement Ar p perpendic- 
ular to the LOS can be written as 



PAr p (Vkick, Z) =P{Ar p >Ar tj L>L t } (Vkick, z) 

xP v kick (ai,a 2 ,g), 



(19) 



where P{Ar p >Ar t ,L>L t } is the probability that the hole is 
shining with resulting luminosity L > L t and is displaced 
by a projected distance larger than the threshold Art. In 
terms of the relevant timescales, 



Pat p (Vkick, z)-- 



t{Ar p >r t ,L>L t } (Vkick, z) 



t 



return 



(20) 



xP v kick (ai,a 2 ,g). 



Figure [l4| shows the probability P/\ r that a MBH 
could be ooserved as a spatially-offset AGN during its 
wandering time assuming Case 1. Along both orienta- 
tions considered here (on and off the plane of the gaseous 
disk), the probability peaks at ~ 10 -6 for optical de- 
tection, Vkick ~ 500 kms -1 , and is higher for the 1:10 
merger where the recoiling hole reaches larger apocenter 
distances. The probability that such a source could be 
observed by Chandra instead drops by nearly an order 
of magnitude. For the 1:1 merger remnant the peak of 
the distribution is ~ 10 -6 in the = 90° case. The dis- 
tribution broadens in the = 0° case peaks at a lower 
value since the MBH spends a considerable fraction of its 
return time in a non-active phase. We have carried out 
the same calculation in Case 2, and find that the likely- 
hood of observing spatially-offset AGNs becomes negligi- 
ble. There are several reasons for this: (1) when spins are 
aligned, the maximum recoil velocity is ~ 200 kms -1 ; 



(2) at low recoil speeds the mass carried by the disk is 
Miisk = M., and the effective recoil velocity is Vkick /2; 
and (3) recoil velocities in the direction perpendicular to 
the gaseous disk are not allowed since = 0. 

7. SUMMARY AND DISCUSSION 

We have studied the orbital evolution of recoiling 
MBHs in post-merger environments via high-resolution 
A^-body + SPH simulations and a semi-analytical model 
that combines the results of merger simulations with a 
prescription for MBH accretion. Our results can be sum- 
marized as follows: 

1. We find little evidence for an observable imprint of 
a recoiling, non-active MBH in the gaseous turbu- 
lent medium of our major merger remnant. In our 
simulation, the total energy deposited by the mov- 
ing hole is only 2% of the change of internal energy 
of the gas, even at high Mach numbers. The contri- 
bution to the thermodynamics of the ambient gas 
becomes significant in minor merger remnants in- 
stead, potentially allowing for an electromagnetic 
signature of MBH recoil. 

2. We have built a library of MBH orbits in three gas 
rich galaxy mergers, allowing the holes to carry a 
punctured accretion disk and the associated BLR 
upon ejection, and to accrete from the ISM when 
the disk is exhausted. In the major merger rem- 
nant, the combination of both deeper central po- 
tential well and drag from high-density gas con- 
fines even MBHs with kick velocities as high as 
1,200 kms -1 to within 1 kpc from the host's cen- 
ter. In minor merger remnants gas drag is weaker, 
the potential is shallower, and smaller recoil speeds 
can displace the hole to larger distances. 

3. We have dissected the orbits of recoiling holes 
into regions where they may be detectable as 
kinematically and/or or spatially offset AGNs. 
Kinematically-offset AGNs carry their accretion 
disk and associated BLR for a timescale M^^/M 
and can be observable as such when their broad 
emission lines have a velocity shift > 600 kms -1 
relative to the narrow-line gas left behind. In the 
major merger remnant, moderate recoil velocities 
(700 — 1,000 kms -1 ) provide the most likely sce- 
narios for the detection of kinematically offset nu- 
clei, as the initial disk can fuel the MBH during 
several pericenter (high- velocity) passages. The 
convolution of the total fractional time spent as a 
kinematically-offset AGN with the probability that 
a coalescing binary receives a given kick velocity 
yields the overall probability of observing recoil- 
ing MBHs as kinematically-offset AGNs. This is 
small, of order ~ 10 -T , even in the optimistic case 
(our "Case 1" ) where the spins of the binary holes 
are randomly oriented relative to the orbital angu- 
lar momentum vector and the orbital plane of the 
binary is randomly oriented relative to the galac- 
tic disk. It is totally negligible when the spin s are 
aligned instead (e.g. jBogdanovic et al. 2007), and 
the maximum recoil speed drops to ~ 200 kms -1 . 
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4. At low redshift, spatially-offset AGN can be ob- 
served with both Chandra and HST/JWST for a 
few Myr in major mergers for high recoil speeds (> 
1,000 kms -1 ). In minor mergers, the detect abil- 
ity window shifts to lower recoil speeds (Vkick < 
500 kms -1 ) and extends to up to a few 100 Myr 
in some cases. Under the Case 1 optimistic sce- 
nario, the probability of observing spatially-offset 
AGNs is larger (by two dex) in the 1:10 merger 
than in the 1:1 merger, since in the former larger 
offsets occur for lower, more likely kicks. The peak 
of the probability distribution (10 -6 ) occurs for 
Vkick — 500 kms -1 and optical detection, and is 
lower by about one dex in the case of X-ray detec- 
tion. 

Our semi-analytical model has been calibrated using a 
limited set of numerical simulations. The choice of a gas 
fraction of 10% for the major merger is rather conserva- 
tive; a larger gas content, quite likely at high redshift, 
would cause increased drag on the MBH and reinforce 
our conclusions. For the 1:10 minor merger, representa- 
tive of some of the most common types of mergers during 
the hierarchical assembly of galaxies, we have assumed 
a larger gas fraction, 30%, as in this case the pairing 
of the two holes and their coalescence is significantly 
faster. Such a choice is still on the conservative side; 
the typical gas fraction in galaxies with masses in the 
range 2 x 10 10 — 2 x lO n M0 (corresponding to the virial 
masses of the primary and secondary galaxies in the 1:10 
merg ers, respectively) is in the range 30 -50% at low red- 
shift flSimon fc Geha||2007| |Leroy et al.||2008| ) and would 
be even higher at early epochs. As for the mass ratio, 
our simulations sample quite well the regime in which 
the dynamical friction timescale is short enough for the 
galaxies to merge in a few orbits. At larger mass ratios 
the satellite galaxies will suffer severe tidal mass losses as 
dynamical friction reduces their pericentre d istance: this 
will increase substantiall y their decay time ( |Colpi et aL 
1999; |Taffoni et aT]|2003[ ) and inhibit the formation of a 
MBH binary. 

A limitation of the present study is that our simula- 
tions do not capture all the phases of the ISM: gas is 
either warm, with temperatures close to 10 4 K, or hot 
(T ~ 10 5 — 10 6 K) as a result of heating by supernovae 
explosions or shocks in the final collision of the merging 
galaxy cores. The cold molecular phase that is believed 
to represent a s ignificant fraction of the m ass of the nu- 
clear disk (e.g. Downes & Solomon 1998) is missing in 
our current calculations. At temperatur es of hundreds 
of ke lvins in the presence of a starburst ( jKlessen et al. 
2007), molecular material may be arranged in a th in and 
dens e sheet the m idplane of the nuclear disk (e.g. |Wada 



[fc Norman 2001). MBHs recoiling along the plane may 
then sutler more dynamical friction, and one may expect 
a greater difference between the fate of in-plane and off- 
plane recoils than seen in our study. 

The finite resolution of our numerical simulations may 
also affect our results. In particular, the resolution of 
our 1:4 merger (200 pc) results in a sharp decline of the 
gas distribution in the remnant nucleus. To quantify this 
effect, we have rescaled the central gas and stellar densi- 
ties to their corresponding values at 200 pc. We find that 
the recoiling MBHs in the 1:4 merger reach smaller apoc- 



enter distances and wander for even shorter timescales. 
The effect is more dramatic for kicks below 200 kms -1 . 
Since MBHs with these velocities are not observable as 
neither kinematically nor spatially-offset AGNs, our de- 
tection probabilities remain largely unaffected. 

Another important limitation of our semi-analytical 
modelling is the neglect of AGN feedback, which acts 
to reduce the efficiency of frictional drag by heating the 
gas surrounding t he recoiling MBH. Recent work by Si- 
jacki et ah] ([2010) shows that in simulations where AGN 
feedback is fully included, the MBH creates a hot low- 
density (instead of a high-density) wake that expands 
away from the hole. The MBH describes a precessing el- 
liptical orbit that tend to circularize, and decays towards 
the center of the remnant on a longer return timescale. 
To get a simple estimate of the effect of AGN feedback 
on our calculations we have assumed in our semi- analytic 
model that at each time- step a fraction equal to 5% of the 
radiative energy produced by the accreting MBH goes di- 
rectly towards increasing the local thermal energy of the 
gas. We find that in dense regions where the hole is 
accreting at high Bondi rates, the reduction of the effec- 
tive Mach number from such energy deposition lengthen 
the return timescale by 15% in the major merger for 
Vkick = 1000 kms -1 and 20% for Vkick = 500 kms -1 in 
the 1:10 merger. The time during which the MBHs are 
observable as spatially offset AGN increases by 6% and 
9% respectively. The overall detection probability how- 
ever, remains low. 

Furthermore, in our detect ability calculations we have 
not included obscuration by dust. Because theoptical 
and X-ray reddening does not depend on the orientation 
of the host galaxy (Winter et al. 2010), it is believed 



that most of the obscuration occurs m a region close to 
the MBH, with optical light being more stro ngly affected 
that X-rays and infrared wavelength s (e.g. |Goulding fc 
Alexanderj|2009 [Winter et al. l| 2010|) . In particular, the 
dToF Koss et al.l 020 



Koss et al.| ^OTPf shows that the Swift 



recent study of 

BAT X-ray survey has revealed roughly 5-6 times more 
dual AGNs than the Sloan Digital Sky Survey (SDSS) 
due to optical extinction and dilution by star formation. 
As a result, the detection probability is further reduced 
for HST/JWST sources. 

While we have pointed out the low likelihood that 
moderate-to-large recoils can be observed, one possibly 
observable consequence of this phenomenon is the ex- 
istence of a population of hot-dust-poor AGNs due to 
the removal of the MBH and associated BLR from the 
surrounding hot dust region. Nearly 10% of spectro- 
scopicaHy confirmed type 1 AGNs samples ( Elvis et al. 
|1994[|Richards et al.|2006||Hao et al.|2010a[ ), and XMM- 
COSMOS (Elvis et al. 2010 in prep) show a relatively 
weak infrared bump, associated with dust emission. The 
number of these objects has been shown to increase with 
redshift, from 6% at z < 2 to 20% at 2 < z < 3.5 ( |Hao 
et al.|2010b[ ). They find that the covering; factor of these 



AGN is only 2-30%, significantly smaller than the 75% 
predicted by the unified model. Since these AGN are in 
the redshift range < z < 4, these hot-dust-free AGN 
are not recent ly born AGN wh ich have not time to form 



a dusty torus ( | Jiang et al.|2010 ). Based on their spectral 
resolution (2 A), a z = 2 source could have a relative 
velocity of 50 — 70 kms -1 . 
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